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1 INTRODUCTION 



ABSTRACT 

We constrain the amplitude of primordial non-Gaussianity in the CMB data taking 
into account the presence of foreground residuals in the maps. We generalise the 
needlet bispectrum estimator marginalizing over the amplitudes of thermal dust, free- 
free and synchrotron templates. We apply our procedure to WMAP 5 year data, 
finding /nl = 38 ± 47 (1 a), while the analysis without marginalization provides 
/nl = 35 ± 42. Splitting the marginalization over each foreground separately, we 
found that the estimates of /nl are positively cross correlated of 17%, 12% with the 
dust and synchrotron respectively, while a negative cross correlation of about —10% 
is found for the free-free component. 

Linde fc Mukhanovl 1200 



There has been considerable activity recently in constraining 
the amount of non-Gaussianity present in CMB data. This 
is principally motivated by theoretical interest into devia- 
tions of primordial fluctuations from Gaussian statistics — 
a natural outcome of several implementations of the infla- 
tionary scenario, which can be used as a tool to rule out 
specific models. More mundanely, non-Gaussianity can also 
be produced by undetected systematics, which may suggest 
problems in the dataset. It could also indicate that elements 
are missing in the standard cosmological model; for exam- 
ple, anomalies which have recently been detected on large 
scales have been put forward as evidence of cosmic defects 
or an anisotropic universe. 

The primordial non-Gaussianity in inflationary mod- 
els is usually characterised t hrough the in troductio n of th e 
/nl parameter (see e. g. iLuol (fl99ll); Heavens! (Il998h ; 
ISpergel fc Goldberg! (|l999l ); iKomatsu fc Spergell (|200lD ) 
which sets the amplitude of the non-linear contribution with 
respect to the leading order: 



20021; iLinde fc Mukhanovl 12006); lAlabidi fc Lvthl 
Mizuno et al.ll2008l;lKhourvll2002l;fsteinhardt fc Turokl 



2006 



2002 



Lehners fc Stein hardt 2008) i t can take much larger val- 



ues (see e.g. iBartolo et al.l l|2004 l for a review). Since 
the release of the first yea r WMAP data (|Bennett et al.l 
120031 ; IKomatsu et~ai] 120031 1 . there has been a drastic re- 
duction of the upper limits of /nl. The most recent con- 
straints coming from WMAP data using d i fferent tech- 



niques c an be found in Komatsu et alj (l2009f) . Smith et al.l 



d2009lh lYadav fc WandeltJ d2008l) ICurto et all d2009l). 



ICurto et all (120091). iPietrobon et ail (120091). IRudiord et ail 
d2009l). IRudiord et all (|2009l ). ISmidt et all (|2009l ) and 
IVielva fc Sand ( |2009l ). Resu lts on /nl with suborbital ex- 
perim ents can be found in dNatoli et al.l 120091; ICurto et al.l 
200Sl). Interestingly, some of these (|Yadav fc Wandeltll2008l ; 



Rudiord et al.ll2~009l ) reported a non-null detection of /nl at 



$(x) = $ L (x) + / NL ($l(x) 2 - ($l(x) 2 » 



(1) 



where < I > (x) is the primordial gravitational potential and 
$l(x) is a Gaussian field. In the following we will fo- 
cus on primordial non-Gaussianity of local type described 
by Eq. [T] Different inflationary scenarios predict differ- 
ent values of /nl; while in the standard slow roll in- 



more than 2a level in the WMAP data, leading to a debate 
on the significance of the signal and whether it may have 
a cosmological origin or be due to residual fo reground or 
instrumental contamination (|Smith et al.ll200"g| ) . 

Non-Gaussianity is also expected from unremoved 
contamination from astrophysical sources, or foregrounds. 
When component separation te ch nique s are applied to 
CMB data (e.g. iMaino et al l d2002l); iTegmark et all 



il2003lh 



Leach et al.l (|200Sf ); 



flationary scenario dGuthl Il98ll ; ISatol Il98ll ; iLindel 1 19821 : 
lAlbrecht fc SteinhardtHl982l), /nl it i s predicted to be of the 



order of unity (see also ( Maldacenal 120031 ; lAcouaviva et al 
120031 )). in several alternative models ( Lvth fc Wands 



Bonaldi et all (120071) ; 
iBernui fc Reboucasl ((20091 )), residual foregrounds re- 
main, and while they are subdominant, they can still be a 
source of non-Gaussianity, and this could be confused with 
a primordial signature and affect the constraints on /nl. 

In th is paper, which is co mplementary to our previ- 
ous work dPietrobon et al.ll2009l ). we aim at generalising the 
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needlet bispectrum estimator to the case in which such fore- 
ground residuals are present in the data. Another marginal- 
ization technique, using the fas t cubic estimator o n the bis- 
pectrum data, can be found in ISmith et all (|2009T ). 

The paper is organised as follows. In section [2] we give a 
brief introduction of the needlets and their bispectrum; sec- 
tion[3]addresses the generalisation of the needlet bispectrum 
estimator in presence of foreground residuals, in section [4] 
we describe the CMB and the foregrounds dataset used and 
show our results, in Section [5] we comment on our findings. 



2 NEEDLET FORMALISM AND GENERAL 
APPLICATIONS 

Here we give a brief summary of the needle t formalism, but 
one should refer to iMarinucci et al.l (|2008t ) and references 
therein for a more com plete discussion. The mathematical 
appro a ch is described in Narcowich et al.l (1200 6 ) , Baldi "et al.l 
l|2006h . lBaldi et all l|2007l ) and lBaldi et al] ^OOa T" 

Needlets are filter functions with some appealing prop- 
erties; in particular they are well localised both in real and 
harmonic space, and are straight forward to implement in 
practice. The analytical function describing a needlet is: 



lM7) = V^EKs^) £ Y tm {*f)Y lm {t jk ). (2) 

l m=—i 

Here, 7 is the direction of a given point on the sphere, {£jk} 
are the cubature points, weighted by Xjk, for a given spatial 
frequency j and location k and b(-) is a filter function de- 
fined in harmonic space. B is a free parameter which, once 
fixed, determines the width of &(•) and determines the range 
of multipoles which appears in Eq. [2] Its value is set ac- 
cording to the angul ar scales one is interest ed in exploring. 
It can be shown (see IMarinucci et al.l (|2008l )) that assuming 
certain properties of the function &(•) (e.g. finite support in £- 
space and smoothness) , the needlet s are quasi-exponentiall y 
localised around the 7 direction (|Narcowich et all l2006h . 
These properties allow needlets to be used as a scale de- 
pendent analyser in both the harmonic and pixel domains. 
For the reader interested in the details of a pr actical imple- 
ment ation of need lets using Healpix pack agqj (|G6rski et al.l 
120051 ). we refer to lPietrobon et aT, I <|2006h . From the opera- 
tional point of view, the quantities of interest are the coef- 
ficients of the expansion on the needlets basis defined as: 



03k 



T(7)^ fc (7)dfi 



J X i k ^ b [^j) E Oim^mfefc)- (3) 
I m=-£ 

These coefficients can be used for reconstructing a map as: 
r(7) = £>klM7) (4) 

j,k 

As with other varieties of wavelets, needlets allow a multi- 
scale analysis, which is usefu l for many kinds of data 
analyses, such as denoising (ISanz et al. (1 19991)) , poin t 
source extraction (see for instance. ICavon et al] l)2000h : 

1 http://healpix.jpl.nasa.gov 



Vielva et al.l d2003l)), asymmetry studies dCruz et al. 2007: 


Pietrobon et all 200S 


: McEwen et al.l 120081: Vielva et all 


2007t iMcEwen et alJ 


2008|) and for cross correlating dif- 



ferent datasets, as for example CMB and large scale sur- 
veys for the detection of the Integrat ed Sachs- Wolfe effect , 
which sets constra ints on dark energy (jPietrobon et al hood : 
I Vielva et all200d ). From the coefficients of Eq.0 we can es- 
timate the binned CMB power spectrum: 



E & 



2 2^+1 

J 47T 



(5) 



(see iFav et al.1 |200Sl ) and IPietrobon et~ai1 (|200d . l2008h for 
an application to real data). It is straightforward to gener- 
alise this for malism to higher order s t atistics like the needlet 
bispectrum dLan fc Marinuccil 120081: IPietrobon et al1l2009l: 
Rudjord et al.l 120091 : IPietrobon et al.l 120091 : IRudiord et al l 
20091 ): 



S 



2J3 ~ AT E 
JV P k 



0j 1 k/3j 2 k/3j 3 k 



J Jl w J2 u 33 



(6) 



Here k spans over the pixels of the region of interest, N p 
is the total number of these pixels and <jj is the stan- 
dard deviation of the needlet coefficients. Eq. [6] has been 
applied to constrain /nl with a needlet analysis, since it 
has b een proven to be proporti onal to the reduced bispec- 
trum (K omatsu fc Spergelll200ll) , and t hus to the non- linear 
parameter f NT , jPietrobon et al]|2009l; IRudiord et al]|2009l ; 
IPietrobon et al-lhoOSl IRudiord et al.ll2009l ): 



S313233 « /nl y^S0fik0?2k0%k + perms) 



(7) 



The needlets bispectrum has several useful properties. It can 
be shown that, given the localisation properties of needlets, 
their bispectrum is not heavily affected by the correla- 
tion introduced by the masking procedure. Another advan- 
tage is that there is no need to calculate Wigner coeffi- 
cients (as a pure bispectrum analysis requires), which is a 
time consuming step, especially when high multipoles are 
considered. These properties translate into a good detec- 
tion power in constraining primordial non-Gaussianity (see 
JPietrobon et all 120091 : IRudiord et all 120091 : IRudiord et al.l 
120091)) comparable with th e fast cubic estim ator introduced 
by Komatsu et all (|2005|) and app lied in iKomatsu et al.l 
(|2009l ) and lYadav fc WandehJ (|2008l ). 



3 METHOD 

In this section we g eneralise the /nl esti mator from the 
needlet bispectrum (|Pietrobon et al.l 120091) in presence of 
residual foregrounds. 

When constraining primordial non-Gaussianity, one can 
estimate /nl by minimising the chi-square: 



X (/nl) 



Y T C' 1 Y, 



(8) 



where Y = Y° bs — (K(/nl)) is the difference between an or- 
dered array of data (Y obB ) and the corresponding theoretical 
prediction Y^/nl). The data can take many forms, such as, 
the values of Minko wski functionals for different thresholds 
(|Hikage et al.l l200d ) , the densities of valleys or hills when 
using the local curvature (|Cabella et al.ll2005|) or. as in this 
case, the values of the needlet bispectrum ( Pietrobon et al.l 
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120091 ; iRudiord et al]|2009l ; iRudiord et ai1l2009l ) . In the pres- 
ence of the expected weak non-Gaussianity, the covariance 
matrix C can be calculated via Gaussian simulations with 

the same power spectr um as the observed data. 

It has been shown (|Pietrobon et al.ll2009l ; IRudiord et al.l 

l2009h that an unbiased estimator for /nl is given by 

EQobs^-i — 1 oth 

L V 5 th C _1 5 th ' 

where n runs over the triplets {,71.72.73} and S^} 1 represents 
the ensemble average of primordial non- Gaussian realisa- 
tions (/nl = 1). Although the process of foreground reduc- 
tion could make things more complicated, we can make the 
minimal assumption that the final map, contaminated by 
foreground residuals, can be modeled as: 

T sim (7) = T G (7) + /nlT ng (7)+7V(7) 

+ a v D(^)+a F F(^) + a s S(j), (10) 

where D, F and S are the thermal dust, free-free emission 
and synchrotron radiation maps, respectively and the noise 
map TV is assumed to be Gaussian. Calculating the needlets 
coefficients from Eq. [10] we obtain 

Pjk = P% + /3j k + /NL/3j*P + aoPfk + a-vPjk + asPjk- (11) 

In Figure [T] we show the needlet coefficients in the case of 
B = 2 and j = 6 for the three foregrounds templates - dust, 
free-free and synchrotron - once they have been masked. 
The maps were converted in thermodynamic temperature 
for each channel and combined to form one single map. Each 
template was divided by a factor 10, which is the level of 
residuals expected. This factor was estimated in harmonic 
space through a Monte Carlo Markov Chain algorithm. The 
foreground reduced map used in this analysis has been pro- 
duced by optimally combining the W, V and Q WMAP 
bands as described below (Eq. 1171 ) 

The angular power spectrum of this optimal CMB map 
has been compared to a linear combination of synchrotron, 
free-free, Galactic dust, CMB and noise. Synchrotron, free- 
free and Galactic dust maps are estimated from the data 
using exactly the same combination in Eq. 1171 while CMB 
and noise are estimated through simulations for each band 
and then combined in the same way as data. 100 realisations 
of CMB maps have then been generated in order to take into 
account the uncertainty due to the cosmic variance. 100 real- 
isations of noise have been created with variance consistent 
with the data and then fitted the amplitude coefficients p 7 
of each component 7, such that 

C j°t = C ? MB + J2p-yCj + Ne, (12) 

7 

where Nt is th e detector noise (a sim ilar procedure has 
been applied in IVeneziani et al.l |2009|) ) . We find p sync h < 
0.1 x 10 _1 , pdust < 1.0 x 10" 1 and p ff < 0.6 x 10 _1 . The 
coefficients p 7 provide an indication of the power spectrum 
contamination percentage by the foreground residuals. The 
effect on the bispectrum is estimated to be of the order of 
~ 10 -3 . 

We can compute the needlet bispectrum expected from 
the model of the model of residuals described in Eq. 1101 by 
substituting Eq. [TT]into the definition of the bispectrum. In 
principle, this will include many terms including triple prod- 
ucts of the various contributions. On average, the product of 




Figure 1. The needlets coefficients with B = 2 and j = 6, which 
are sensitive to 64 < i < 127, or angular scales of a few degrees. 
The maps shown are for thermal dust (top), free- free, synchrotron 
and the assumed primordial non-Gaussian signal for /nl = 50 
(bottom). The Kq75 mask has been applied. 

the Gaussian contributions should be zero, and the lowest 
order contribution from the primordial non-Gaussianity is 
that shown in Eq. [7] We also expect on average that cross 
terms like /3 NG /3 :i /3 l2 (I={D,F,S} will also be zero, as the 
foregrounds should be uncorrelated with the primordial sig- 
nal. 

However, cross terms which include different fore- 
grounds, of the form (3 11 f3 l2 /3 l3 , could potentially be non- 
zero, particularly if there are physical reasons to expect cor- 
relations. (For example, both synchrotron and free-free are 
sensitive to the free electron density.) Even if such terms 
are zero on average, they will not be zero for specific reali- 
sations of the foregrounds; indeed if we think we know the 
foreground templates, we can calculate these terms directly. 
For simplicity here, we ignore such terms and focus only on 
the 'auto-bispectra.' This has the advantage that solving for 
the best fit amplitudes involves solving a simple set of lin- 
ear equations. (Cross terms make the equations non-linear, 
requiring a numerical solution.) 

With this assumption, the theoretical needlet bispec- 
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Figure 2. Bispcctrum amplitudes for the primordial non- 
Gaussian signal (/nl = 50), and the three WMAP data fore- 
ground templates suppressed by a factor 1000, roughly corre- 
sponding to the expected signal of the residual foregrounds. 



trum (5) in presence of foreground residuals can be written 



3233 ' 

(13) 



where 



S 



k 



aG oG oNG 
Pj 1 kPj 2 kPj 3 k 



' 3\ w 32 u 33 



Pjl kPj 2 kPj 3 k 



J 31^ 32^33 



+ perms; 



(14) 



For ease of notation, we have denned the foreground ampli- 
tudes as Ai = a? . 

In Figure [5] we show the bispectrum for the primor- 
dial non-Gaussianity with /nl = 50 and for the foreground 
templates of WMAP data normalised according to the pre- 
scription discussed previously. The bispectrum is shown as a 
function of the quantity X — + l/jija + l/j3 32, which 

orders them roughly according to their variance (as can be 
seen in Figure [3]) 

The x 2 expression (Eq. [S} can be easily generalised as 
follows: 



x 2 = {s; bs - sm^) - sUAi))c^, 
(^v s -s m ,(/ nl )-s;,(ao) 



T 



(15) 



where i refers to the i-th foreground template with the Ein- 
stein summation convention assumed, and is the needlet 
bispectrum of the given foreground. Minimising Eq. ll5l with 
respect to x = (/nl, An, Af,As) we obtain: 



/nl 



u y 3 



oNGf-l qNG 



.Z_^ oNG^-lcNG K 



Cl/^f— 1 Qobs 

= „,„-i g i - 2^ 



ol/^— 1 qK 

— ; 1 — r-A* 



(16) 



where in each case the sum implicitly goes over the other 
signals K = {NG,D,F,S}. 

The solution of the previous system provides us with the 
estimates of /nl with the needlet bispectrum in presence of 



foreground residuals. In the following we present data and 
simulations where this estimator has been applied. 



4 DATA SET, SIMULATIONS AND RESULTS 

In the following, the needlets of the simulations and data will 
be calculated for B=2 unless specified othe rwise. We used 
the p ublicly available WMAP 5-year dat£0 | Hins haw et al.l 
I2009T ). The CMB map has been obtained combining the fore- 
ground reduced maps for the channel Q,V,W according to: 



T obs ( 7 ) 



E np bB (7U 



(17) 



where n b s (7) is the number of hits for a given point on 
the sphere (7), a c h is the nominal sensitivity for the rel- 
ative channel ch and the superscript red denotes the re- 
duced dataset. We considered the WMAP data templates 
of dust, free- free and synchrotron emission at the resolu- 
tion of Af s ide = 256 corrected by the conversion factors 
from antenna to thermodynamic temperature of the respec- 
tive channel; we then combined them as in Eq. [T7] to have 
one single map for each template. The data and foreground 
templates were then masked with Kq75 mask (downgraded 
to the same resolution as well) which covers the dominant 
Galactic emission over roughly the 30% of the sky. We ob- 
tained our final constraints on /nl by applying the estimator 
in its improved fashion to the WMAP 5-year data where: 

• the covariance matrix C _1 was calibrated over 20,000 
Gaussian simulations; 

• the needlet foreground bispectra 5d , Sf , Ss were calcu- 
lated on the templates described above; 

• the primordial needlet bispectrum 5* NG was calculated 
using over 100 primordial non-Gaussian maps (|Liguori et al.l 
l2007h convolved with the beams and combined as done for 
the Gaussian simulations (see also lElsner &i Wandeltl (|2009h 
for an alternative algorithm to generate temperature and 
polarization primordial non Gaussian maps). 

Figure [3] shows the data we used together with the av- 
erage and the standard deviation derived from simulations. 
The error bars on /nl were computed through the distri- 
bution of its estimates for 20,000 Gaussian realizations. In 
Figure [4] we show our results, where in each panel we show 
the estimate of /nl with and without marginalizing over 
the foreground templates. Without marginalizing, we found 
/nl = 38 ± 42, 85 at la and 2a, respectively. The marginal- 
ization over all the foregrounds brings the constraint to 
/nl = 41 ±47, 95 with an increase of the error bars of about 
10% and a positive shift of the mean value. Although the 
larger error bars mean that we should not attach too much 
significance to the higher /nl value, it would be interesting 
to investigate whether this could partiall y explain some pos- 
itive d e tection present i n liter ature (e. g. lYadav fc Wandeltl 
(|2008l ): iRudiord et all lj2009h ). Further information can be 
obtained by looking at the estimates of /nl marginalizing 
over the three foregrounds separately. The enlargement of 
the error bars due to each foreground is of the same order of 
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Figure 4. One dimensional likelihoods for the /nl parameter. Solid lines represent marginalized likelihoods, whereas the dashed ones 
refer to the slice j4f oreg = 0. The top left panel shows the analysis including all the three foregrounds, while the others are derived when 
a single foreground — dust (top right), free- free (bottom left) and synchrotron (bottom right) emission - is allowed to vary. 



ICO - 



2C0 



-2C0 - 



-100 




Figure 3. Needlets bispectrum measured from the WMAP 5- 
year temperature map. The shaded region marks the 68% confi- 
dence level computed from 20,000 Gaussian realisations. 



magnitude, but the shift of the estimate seems to be mostly 
due to the dust component. 

As a further check we carried out a Fishe r analy sis. (For 
a quick-start guide to Fisher matrix see lCoej l|2009l ).) So far 
we applied a Monte Carlo approach to find a good estimate 
of /nl in presence of foreground contamination and quantify 
the scatter around this evaluation. 



If the likelihood is close to Gaussian, the scatter can be 
well approximated by the Fisher approach, as the Cramer- 
Rao bound ag 1/(.F _1 )V 2 becomes close to an equality. 
The Fisher matrix is defined as: 



da 



-c~ 



db 



(18) 



where S is the signal described in Eg. 1131 and p, p , a and 
b run over the triplets {ji , J2 j J3 } an d the parameter set 
{/nl, An, Af, As}. In detail we have: 



F 



S NG c — l S NG 
S U C -1 S NG 
s P c _l gJ VG 



S NG C 
S D C~ 
S F C" 

s s c~ 



S NG C - 

s D c- 
s s c~ 



S NG 
S D C 
S F C 

s s c 



c~ 



where the sum over the triplets is implicit. 

The marginalized error on /nl arising from this Fisher 
analysis is A/nl = 45, very close to the limits obtained 
with the Monte Carlo approach; this confirms the efficiency 
of our estimation method. In Fig. [5] we show the ler and 2a 
Fisher error ellipses together with the Markov chain output; 
here each plot presents the significance region in the case 
where the other parameters are fixed at their fiducial value. 
It can be seen that the free-free component seems slightly 
anti-correlated with /nl, the cross-correlation coefficient 
being p = ^S NG C- 1 S F /(S NG C- 1 S NG x S F C^S F )^ 2 = 
—0.10. The other foreground components, i.e. dust and syn- 
chrotron, show a positive correlation, demonstrated by the 
orientation of the ellipses, with p = 0.17 and p = 0.12, re- 
spectively. We can also calculate the cross-correlation be- 
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tween the various foreground amplitudes based on their 
non-Gaussianity: the dust-synchrotron correlation is p — 
—0.87, the dust-free-free correlation is p — —0.29 and the 
synchrotron-free-free correlation is p — —0.58. The scatter 
of the over-plotted distribution of Gaussian simulations on 
the same Ai-/nl plane is in excellent agreement with the 
ellipses, which confirms again the consistency of our proce- 
dure. We can also use the scatter in the foreground ampli- 
tudes to estimate limits on the foreground contamination 
using the bispectra; these are generally the same order as 
the limits arising from the fitting of the power spectrum 
discussed above. 

As a check of our assumption that we are not strongly 
affected by potential cross-correlations among the fore- 
grounds templates, we also considered the three foreground 
signals combined in a single foreground template. This ef- 
fectively maps to a perfect covariance between the various 
foreground amplitudes. The analysis for this combination 
template is identical to the individual foreground analyses 
discussed above. We use A s to describe the amplitude of 
the the needlet bispectrum of the total foreground map, ob- 
tained summing the three templates. We found /nl = 37±43 
at 1 a confidence level, very close to the result when the fore- 
grounds are included individually and assumed to be uncor- 
rected. We also plot the two-dimension scatter plot in the 
plane fyh-A s in Fig. [5] The correlation results p ~ 0.14, 
in good agreement with the sum of the three correlations 
of thermal dust, free-free and synchrotron emission, con- 
forming the reliability of the approach we followed. Fi- 
nally, we also applied our improved estimator to the raw 
WMAP 5-year temperature maps, where no foreground re- 
moval was attempted. In this case, the estimate of the pa- 
rameter /nl decreases to 20: this is consistent with the over- 
all positive value we obtain for the correlation between fore- 
grounds and primordial non-Ga ussianity and is also co nsis- 
tent with earlier work, s uch as iKomatsu et~aH (|2009h and 
lYadav fc WandeltJ (|2008T ). 



5 CONCLUSIONS 

In this paper we have presented a procedure to marginalize 
the residual foregrounds when estimating /nl in the needlet 
bispectrum framework. However it is important to stress 
that this algorithm does not strictly rely on needlets prop- 
erties and it can be easily applied to any linear estimator. 
With the foreground marginalization, we found, for WMAP 
5-year data, that the error bars are enlarged by about 10% 
with respect to the estimate obtained without marginaliz- 
ing. Foreground residuals can have different effects when dif- 
ferent estimators are used to characterised primordial non- 
Gaussianity. Comparing with other foregr ound analyses, our 
result s seem to go in the direction of lYadav fc Wandeltl 
(|2008l ). where they arg ued that foregroun ds negatively bi- 
ased the /nl ■ However, ISmith et al. | (j2009l ) draw other con- 
clusions showing that the sign of the biasing could depend 
on the choice of weighting in the method of estimating /nl- 
Since needlets coefficients are essentially a rebinning of the 
filtered harmonic coefficients ai m , this could explain the dif- 
ferent behaviours. This issue probably needs to be exam- 
ined separately for each test of non-Gaussianity, since dif- 
ferent tests act differently on different spaces (e.g. harmonic, 



pixel, wavelet), suggesting that the influence of foregrounds 
on estimating /nl is not unique. 

All this reinforces the argument that a careful analy- 
sis with different tests of non-Gaussianity is crucial to dis- 
criminate between primordial non-Gaussianity and spuri- 
ous effects. Our procedure could be improved further by 
incorporating any covariance which may arise in the fore- 
ground subtraction methods, and by fully considering the 
effects arising from the potential correlations between dif- 
ferent foreground maps. Such analyses will be essential for 
modelling the non-Gaussianity of future experiments like 
Planck, where the error bars on /nl are expected drast i- 
cally reduced (A/nl ~ 3 - 5 IKomatsu fc Spergell (|200ll ); 
iBabich fc Zaldarriagal |200i)), making the uncertainties in- 
troduced by foregrounds extremely relevant. 
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